
Job: betamc.std.prg (/home/gov/faculty/ppaolino/parr)
Sun Jun 17 00:34:23 2001, pid=16176, uid=1141(ppaolino)

cheng alg used if cheng=1, johnk otherwise, cheng=       1.0000000 
parm1 =        1.5000000       0.30000000 
parm2 =      -0.40000000      -0.30000000 
p =        4.7000044 
q =       0.69877941 
true first difference (mean)=      0.13982664 
true first difference (variance)=    -0.023162353 

descriptive stats for standard beta parameters

       1.6546775       0.30402866 
      0.28544900       0.35531738 
     -0.28629594       0.23338207 
     -0.31469110       0.25164987 
 
standard beta mse=     0.018920615 
 
alternative beta mse=     0.018957440 
 
ols mse=     0.019107655 
 
logit trans ols mse=     0.023023084 
 
het ols mse=     0.019906151 
 
logit trans het ols mse=     0.022502188 
 
log trans het ols mse=     0.019998214 
 
first differences for mean and variance for std beta (sim)
      0.13493333      0.054503486 
    -0.019564358      0.014692799 
 
first differences for mean and variance for std beta (nosim)
      0.13521405      0.053907728 
    -0.019759440      0.014205716 
 
first differences for mean and variance for alt beta (sim)
      0.13296937      0.053412215 
    -0.019171512      0.014442103 
 
first differences for mean and variance for alt beta (nosim)
      0.13329373      0.052818661 
    -0.019284862      0.013842664 
 
first differences for mean ols
      0.14328512      0.056680301 
 
first differences for mean ols (logit trans)
      0.12157761      0.051011427 
 
first differences for mean and variance for het normal
      0.11724352      0.043321485 
    -0.025487415      0.020105222 
 
first differences for mean and variance for logit trans.
      0.11524966      0.048336816 
    -0.025972703      0.020230919 
 
first differences for mean ols (log trans)
      0.16850329      0.075504608 
 
efficiency of standard beta vs. alternative beta=       98.367318 
 
efficiency of beta vs. ols=       104.95504 
 
efficiency of beta vs. ols (log trans)=       100.13971 
 
efficiency of standard beta vs. het normal=       90.305276 
 
efficiency of standard beta vs. het normal (log trans)=       100.23440 
 
efficiency of standard beta vs. OLS (log trans)=       149.28781 
 
efficiency of standard beta vs. alternative beta=       98.411788 
 
efficiency of standard beta vs. het normal=       138.54973 
 
efficiency of standard beta vs. het normal (log trans)=       139.82319 
 
standard beta parms' overconfidence=
       103.44700 
       111.10735 
       104.82334 
       103.44209 
 
alternative beta parms' overconfidence=
       104.99869 
       107.46131 
       103.93647 
       109.10919 
 
ols parms' overconfidence=
       97.002034 
       106.61016 
 
ols (log. trans) parms' overconfidence=
       99.283546 
       101.27241 
 
ols (log trans) parms' overconfidence=
       96.946204 
       117.34725 
 
heteroskedastic normal parms' overconfidence=
       107.52406 
       89.155142 
       155.34180 
       159.79983 
 
heteroskedastic log. trans. parms' overconfidence=
       103.82384 
       105.84823 
       141.91355 
       135.75150 
 
standard beta's mean first difference overconfidence
       102.57133 
 
standard beta's variance first difference overconfidence
       96.359026 
 
alternative beta's mean first difference overconfidence
       100.75851 
 
alternative beta's variance first difference overconfidence
       94.347472 
 
check standard beta's percentage within 95% conf int

      0.93700000 
      0.92100000 
      0.94000000 
      0.94500000 
 
check alternative beta's percentage within 95% conf int

      0.93200000 
      0.92300000 
      0.93600000 
      0.92800000 
 
check ols percentage within 95% conf int

      0.93900000 
      0.91900000 
 
check logit trans/ols percentage within 95% conf int

      0.93400000 
      0.93500000 
 
check het ols percentage within 95% conf int

      0.90700000 
      0.96400000 
      0.79000000 
      0.78800000 
 
check logit trans/ols percentage within 95% conf int

      0.92700000 
      0.92400000 
      0.85100000 
      0.83900000 
 
check logit trans/ols percentage within 95% conf int

      0.91000000 
      0.86600000 
 
check percentage of sig parms (standard beta)
      0.75500000 
      0.28800000 
 
check percentage of sig parms (alternate beta)

       1.0000000 
      0.76700000 
       1.0000000 
      0.13300000 
 
check percentage of sig parms (basic ols)

       1.0000000 
      0.77300000 
 
check percentage of sig parms (logit trans ols)

       1.0000000 
      0.79100000 
 
check percentage of sig parms (het ols)

       1.0000000 
      0.72800000 
       1.0000000 
      0.61400000 
 
check percentage of sig parms (het ols - logit trans)

       1.0000000 
      0.79400000 
      0.92900000 
      0.35800000 
 
check percentage of sig parms (ols -- log tranform)

       1.0000000 
      0.75100000 
 
descriptive stats for y's for successful runs

      0.85635109 
      0.15776944 
      0.39671384 
      0.99921075 
descriptive stats for y's for failed runs
       0.0000000 
descriptive stats for x's for successful runs

     0.011698864 
       1.0154969 
      -1.8827054 
       1.7975682 
descriptive stats for x's for failed runs
       0.0000000 
number of failed runs
       0.0000000 

code used to generate the above results:

library cml;
#include cml.ext;
CMLset;

fails=0;
crash=0;
ystatf=0;
xstatf=0;

/* load x matrix based upon number of observations */

x = ones(30,1)~rndn(30,1);
n = rows(x);

/* produce matrices for simulated first differences */

meanx=meanc(x); @ create vector of means for the indep vars @
stdx=stdc(x);   @ create vector of stdev for the indep vars @
minx=minc(x);
maxx=maxc(x);

/* start obtaining the first difference vectors (means) */

sw=zeros(1,cols(x)-1)|eye(cols(x)-1);
xl=meanx-stdx.*sw;           @ these are kXk-1 vectors @
xh=meanx+stdx.*sw;       
xmin=meanx.*(1-sw)+minx.*sw;  
xmax=meanx.*(1-sw)+maxx.*sw;

/* start creating Monte Carlo data */
k = cols(x);
parm1 = {1.5, 0.3};
parm2 = {-0.4, -0.3};
cheng=1;
print "cheng alg used if cheng=1, johnk otherwise, cheng=";;cheng; 
p = exp(x*parm1);
q = exp(x*parm2);

/* determine true first differences */

/* get  alpha and beta values with the j^th independent var 
   at +/- 1 sd and min/max values, while other x's held at their means */
         alt=exp(xl'*parm1);                          @ nxk-1 matrix @
         blt=exp(xl'*parm2);
         aht=exp(xh'*parm1);
         bht=exp(xh'*parm2);

/* get predicted values of mean and variance with the j^th independent var
   at +/- sd and min/max values, while other x's held at their means */
         meanlt=alt./(alt+blt);                    @ nxk-1 matrix @
         meanht=aht./(aht+bht);
         varlt=(alt.*blt)./(((alt+blt)^2).*(alt+blt+1));
         varht=(aht.*bht)./(((aht+bht)^2).*(aht+bht+1));

/* produce first differences for the predicted mean and variance with 
   the j^th indep var at +/- sd, while other x's held at their means */
         fdmt=meanht-meanlt;
         fdvt=varht-varlt;

/* print data parameters for simulation */

print "parm1 = ";; parm1';
print "parm2 = ";; parm2';
print "p = ";; meanc(exp(x*parm1));
print "q = ";; meanc(exp(x*parm2));
print "true first difference (mean)=";; fdmt;
print "true first difference (variance)=";; fdvt;

/* start Monte Carlo trials */
trial=1;
do while trial<=1000; 

print "MC trial=";; trial;

   if trial==1;
      ret=0;
      if fails==0;
         bt=0; 
      endif;
   endif;

   if fails==0;
   success=0;

/* comment out random number generator if p,q<1 */
if cheng==1;
/* create J&K (Chen) alpha */
ralpha=p+q;

/* create J&K (Chen) beta */
p1=p.>1;
q1=q.>1;
rbeta=real((p1.*q1).*(sqrt((ralpha-2)./(2.*p.*q-ralpha))) +
      (1-p1.*q1).*(1/(minc((p'|q')))));

/* create J&K (Chen) gamma */
rgamma=p+1/rbeta;

/* step 1 in J&K (Chen) RNG */
u1=rndu(n,1);
u2=rndu(n,1);
v=rbeta.*ln(u1./(1-u1));
w=p.*exp(v);

/* step 2 in J&K (Chen) RNG */

done=0;
      do while done<n;
         d=(ralpha.*ln(ralpha./(q+w))+rgamma.*v-1.3862944).<ln((u1^2).*u2); 
         u1=(1-d).*u1 + d.*rndu(n,1);
         u2=(1-d).*u2 + d.*rndu(n,1);
         v=(1-d).*v + d.*rbeta.*ln(u1./(1-u1));
         w=(1-d).*w + d.*p.*exp(v);
         done=n-sumc(d);
      endo;

/* step 3 in J&K (Chen) RNG */
y=w./(q+w); 
endif;

/* discretize y 

y1=yc.>0;
y2=yc.>0.14;
y3=yc.>0.28;
y4=yc.>0.42;
y5=yc.>0.56;
y6=yc.>0.70;
y7=yc.>0.84;

y=(y1+y2+y3+y4+y5+y6+y7)/8; */

/* comment out random number generator if p or q>1 */
if cheng/=1;
/* step 1 in Johnk's RNG */

u1=rndu(n,1);
u2=rndu(n,1);

v1=u1^(1/p);
v2=u2^(1/q);

/* step 2 in Johnk's RNG */

w=v1+v2;

/* step 3 in Johnk's RNG */

done=0;
      do while done<n;
         d=w.>1;
         u1=(1-d).*u1 + d.*rndu(n,1);
         u2=(1-d).*u2 + d.*rndu(n,1);
         v1=(1-d).*v1 + d.*u1^(1/p);
         v2=(1-d).*v2 + d.*u2^(1/q);
         w=v1+v2;
         done=n-sumc(d);
      endo;

/* step 4 in Johnk's RNG */
y=v1./w;
endif;

/* cheat function to re-set Ys that equal 1 --- would choke llik function
   otherwise */
d=y.==1;
print "cheat=";;sumc(d);
y=(1-d).*y + d.*.99999;  


/* return the monte carlo data set */
    data = y~x[.,2:cols(x)];

   endif;

   if fails==0;

      do while success==0 and fails<=30;

/* maximum likelihood function for estimating parameters alpha and beta */
__title=" ** A Beta Model, Analytic Solution (Standard Estimation) **";

stval=1.0|ones(k-1,1)*0.5|1.0|ones(k-1,1)*0.4;

/* for failed runs, use starting values of crashed iteration */

       if fails>=1 and fails<=40;
          stval=(bt')/2;
          print "retrying run";
       endif;

print "mean y, max y, min y, std y";
meanc(y)~maxc(y)~minc(y)~stdc(y);

proc psi(x);
  local p;
  x=x+6;
  p=1/(x.*x);
  p=(((0.004166666666667*p-0.003968253986254).*p+
      0.008333333333333).*p-0.083333333333333).*p;
  p=p+ln(x)-0.5/x-1/(x-1)-1/(x-2)-1/(x-3)-1/(x-4)-1/(x-5)-1/(x-6);
retp(p);
endp;

proc gradproc(theta,ds);
    local y,x,z,be,h,xb,p,q,per1,gr1,gr2,grad;
		     
    y=ds[.,1];
    x=ones(rows(y),1)~ds[.,2:cols(ds)]; 
    be=theta[1:cols(x),1];                 @ parameters for alpha @
    h=trimr(theta,rows(be),0);            @ parameters for beta @
    p=exp(x*be);                 @ expression for alpha @
    q=exp(x*h);                  @ expression for beta @

    per1=p+q;
    gr1= psi(p+q).*x.*p-psi(p).*x.*p+x.*p.*ln(y);
    gr2= psi(p+q).*x.*q-psi(q).*x.*q+x.*q.*ln(1-y);
    grad=gr1~gr2;
    retp(grad);
    endp;

proc loglik(theta,ds);
    local y,x,z,be,h,xb,p,q,l1,l2,l3,l4,l5,llik;
		     
    y=ds[.,1];
    x=ones(rows(y),1)~ds[.,2:cols(ds)]; 
    be=theta[1:cols(x),1];                 @ parameters for alpha @
    h=trimr(theta,rows(be),0);            @ parameters for beta @
    p=exp(x*be);                 @ expression for alpha @
    q=exp(x*h);                  @ expression for beta @

    if (sumc((p.>0)))<n;
       llik=miss(0,0); print "run failed"; 
    elseif (sumc((q.>0)))<n;
       llik=miss(0,0); print "run failed"; 
    elseif (sumc((y.>0)))<n;
       llik=miss(0,0); print "run failed"; 
    elseif (sumc(((1-y).>0)))<n;
       llik=miss(0,0); print "run failed"; 
    else;
       llik=lngm(p+q)-lngm(p)-lngm(q)+(p-1).*ln(y)+
         (q-1).*ln(1-y);  @ log-likelihood @
    endif;
    retp(llik);endp;

_cml_GradProc=&gradproc;
_cml_GradCheck=1;
_cml_Options = { none }; 
_cml_MaxIters=400;

   if fails.>=5 and fails.<=10;
     _cml_DirTol=.000001;
   endif;
   if fails.>10 and fails.<=15;
     _cml_DirTol=.0000001;
   endif;
   if fails.>15 and fails.<=20;
     _cml_DirTol=.00000001;
   endif;
   if fails.>20;
     _cml_DirTol=.00001;
   endif;

    {b,logl,g,vc,ret}=CMLprt(CML(data,0,&loglik,stval));

/* check that parameters are reasonable */
   bt=(b');

/* determine success or failure of the trial */
   cor=vc[1,1]/vc[1,1];
         if ret>=1 or cor/=1;
            fails=fails+1;
         elseif ret==0 and cor==1;
            success=1; fails=0;
         endif;
      endo;
   endif;

   if fails==0;
      success=0; bta=0;
      do while success==0 and fails<=30;

CMLset;
/* Alternative beta estimation */
__title=" ** A Beta Model, Numerical Solution (Alternative Estimation) **";

stvala=0.0|ones(k-1,1)|2.0|ones(k-1,1)*0.4;

/* for failed runs, use starting values of crashed iteration */

      if fails>=1 and fails<=30;
         stvala=(bta')/2;
         print "retrying run";
      endif;

proc gradpa(theta,ds);
    local y,one,x,z,be,h,xb,zh,zh1,e,per1,gr1,gr2,grad;
    y=ds[.,1];	     
    x=ones(rows(y),1)~ds[.,2:cols(ds)]; 
    z=x;
    be=theta[1:cols(x),1];
    h=trimr(theta,rows(be),0);      
    xb=exp(x*be);
    zh=exp(z*h);
    zh1=zh+1;
    e=xb./(1+xb);
    per1=(-x.*e+x.*(e^2)).*zh1;

    gr1=psi(e.*zh1+(1-e).*zh1-1).*(x.*e.*zh1-(e^2).*zh1.*x+per1)-
        psi(e.*zh1-e).*(x.*e.*zh1-(e^2).*zh1.*x-x.*e+(e^2).*x)-
        psi((1-e).*zh1-1+e).*(per1+x.*e-(e^2).*x)+
        (x.*e.*zh1-(e^2).*zh1.*x-x.*e+(e^2).*x).*ln(y)+
        (per1+x.*e-(e^2).*x).*ln(1-y);
    gr2=psi(e.*zh1+(1-e).*zh1-1).*(e.*z.*zh+(1-e).*z.*zh)-
        (psi(e.*zh1-e).*xb.*z.*zh)./(1+xb)-
        psi((1-e).*zh1-1+e).*(1-e).*z.*zh+xb.*z.*zh.*ln(y)./(1+xb)+
        (1-e).*z.*zh.*ln(1-y);
    grad=gr1~gr2;
    retp(grad);endp;

proc loglika(theta,ds);
    local y,one,x,z,be,h,xb,zh,e,v,p,q,llik,ds1;
		     
    y=ds[.,1];	     
    x=ones(rows(y),1)~ds[.,2:cols(ds)]; 
    z=x;
    be=theta[1:cols(x),1];
    h=trimr(theta,rows(be),0);      
    xb=exp(x*be);
    zh=exp(z*h);
    e=xb./(1+xb);
    v=e.*(1-e).*(1/(zh+1));
    p=((e^2).*(1-e))./v-e;
    q=(e.*(1-e)^2)./v-(1-e);
    if (sumc((p.>0)))<n;
       llik=miss(0,0); print "run failed"; 
    elseif (sumc((q.>0)))<n;
       llik=miss(0,0); print "run failed"; 
    elseif (sumc((y.>0)))<n;
       llik=miss(0,0); print "run failed"; 
    elseif (sumc(((1-y).>0)))<n;
       llik=miss(0,0); print "run failed"; 
    else;
       llik=lngm(p+q)-lngm(p)-lngm(q)+(p-1).*ln(y)+
         (q-1).*ln(1-y);  @ log-likelihood @
    endif;
    retp(llik);endp;

_cml_GradProc=&gradpa;
_cml_GradCheck=1;
_cml_Options = { none }; 
_cml_MaxIters=400; 
    {ba,logl,g,vca,ret}=CMLprt(CML(data,0,&loglika,stvala));

/* check that parameters are reasonable */
   bta=(ba');

/* determine success or failure of the trial */
   cora=vca[1,1]/vca[1,1];
         if ret>=1 or cora/=1;
            fails=fails+1;
         elseif ret==0 and cora==1;
            success=1; fails=0;
         endif;
      endo;
   endif;

/* linear heteroskedastic estimation */

   if fails==0;
      success=0; bth=0;
      do while success==0 and fails<=30;

CMLset;

__title=" ** A Linear-Normal Model, Numerical Solution **";

stvalh=inv(x'*x)*x'*y|(-0.5*ones(cols(x),1));
       
/* for failed runs, use starting values of crashed iteration */

      if fails>=1 and fails<=10;
         stvalh=(bth')/2;
         print "retrying run";
      endif;  

proc gradph(theta,ds);
    local y,one,x,z,be,h,xb,zh,gr1,gr2,grad;
    y=ds[.,1];	     
    x=ones(rows(y),1)~ds[.,2:cols(ds)]; 
    z=x;
    be=theta[1:cols(x),1];
    h=trimr(theta,rows(be),0);      
    xb=x*be;
    zh=exp(z*h);
    gr1=(y-xb).*x./zh;
    gr2=-0.5*(x-(x.*(y-xb)^2)./zh);
    grad=gr1~gr2;
    retp(grad);endp;

proc loglikh(b,ds);                     
    local xb,x1,y,llik,b1,g1,s,s2,u,ds1,one;
    y=ds[.,1];    /* dependent variable */
    one=ones(rows(ds),1);
    x1=one~ds[.,2:cols(ds)];
    b1=b[1:cols(x1),1];   /* partition off parms for E(y) */
    xb=x1*b1;
    s=trimr(b,rows(b1),0);      /* partition off parms for V(e) */
    g1=one~ds[.,2:cols(ds)];
    s2=exp(g1*s);
    u=y-xb;             
    llik=-0.5*(ln(s2)+u.*u./s2);  /* likelihood function */
    retp(llik);endp;

_cml_GradProc=&gradph;
_cml_GradCheck=1;
_cml_Options = { none };
_cml_DirTol=0.000001; 
_cml_MaxIters=500;

    {blh,logl,g,vch,ret}=CMLprt(CML(data,0,&loglikh,stvalh));

/* check that Monte carlo parameters are reasonable */
bth=(blh');

/* determine success or failure of the trial */
   corh=vch[1,1]/vch[1,1];
       if ret>=1 or corh/=1;
          fails=fails+1;
          elseif ret==0 and corh==1;
          success=1; fails=0;
       endif;
    endo;
endif;

/* linear heteroskedastic estimation (with logit tranformation) */

   if fails==0;
      success=0; btl=0;
      do while success==0 and fails<=30;

CMLset;

__title=" ** A Linear-Normal Model, Numerical Solution **";

stvall=inv(x'*x)*x'*(ln(y./(1-y)))|(-0.5*ones(cols(x),1));
       
/* for failed runs, use starting values of crashed iteration */

      if fails>=1 and fails<=30;
         stvall=(btl')/2;
         print "retrying run";
      endif;  

proc gradpl(theta,ds);
    local y,one,x,z,be,h,xb,zh,gr1,gr2,grad;
    y=ln(ds[.,1]./(1-ds[.,1]));
    x=ones(rows(y),1)~ds[.,2:cols(ds)]; 
    z=x;
    be=theta[1:cols(x),1];
    h=trimr(theta,rows(be),0);      
    xb=x*be;
    zh=exp(z*h);
    gr1=(y-xb).*x./zh;
    gr2=-0.5*(x-(x.*(y-xb)^2)./zh);
    grad=gr1~gr2;
    retp(grad);endp;

proc loglikl(b,ds);                     
    local xb,x1,y,llik,b1,g1,s,s2,u,ds1,one;
    y=ln(ds[.,1]./(1-ds[.,1]));    /* dependent variable */
    one=ones(rows(ds),1);
    x1=one~ds[.,2:cols(ds)];
    b1=b[1:cols(x1),1];   /* partition off parms for E(y) */
    xb=x1*b1;
    s=trimr(b,rows(b1),0);      /* partition off parms for V(e) */
    g1=one~ds[.,2:cols(ds)];
    s2=exp(g1*s);
    u=y-xb;             
    llik=-0.5*(ln(s2)+u.*u./s2);  /* likelihood function */
    retp(llik);endp;

_cml_GradProc=&gradpl;
_cml_GradCheck=1;
_cml_Options = { none };
_cml_DirTol=0.000001; 
_cml_MaxIters=500;

    {bll,logl,g,vcl,ret}=CMLprt(CML(data,0,&loglikl,stvall));

/* check that Monte carlo parameters are reasonable */
btl=(bll');

/* determine success or failure of the trial */
   corl=vcl[1,1]/vcl[1,1];
      if ret>=1 or corl/=1;
         fails=fails+1;
      elseif ret==0 and corl==1;
         success=1; fails=0;
      endif;
   endo;
endif;

/* if trial has failed 30 times */

   if fails==30;
      print "run failed after 30 attempts";
      crash=crash+1;

      if ystatf==0;
         ystatf=(meanc(y)~stdc(y)~minc(y)~maxc(y));
      elseif ystatf/=0;
         ystatf=ystatf|(meanc(y)~stdc(y)~minc(y)~maxc(y));
      endif;

      if xstatf==0;
         xstatf=(meanc(x[.,2])|stdc(x[.,2])|minc(x[.,2])|maxc(x[.,2]));
      elseif xstatf/=0;
         xstatf=xstatf~(meanc(x[.,2])|stdc(x[.,2])|minc(x[.,2])|
                 maxc(x[.,2]));
      endif;
      fails=0;
      clear y;
   endif;

/* collect data for a successful trial */
   if success==1;

/* collect y stats for a successful trial */
      if trial==1;
         ystats=meanc(y)~stdc(y)~minc(y)~maxc(y);
         xstats=(meanc(x[.,2:cols(x)])|stdc(x[.,2:cols(x)])|
                minc(x[.,2:cols(x)])|maxc(x[.,2:cols(x)]));
         bsum=b;
         sesum=sqrt(diag(vc));
         bsuma=ba;
         sesuma=sqrt(diag(vca));
         bsumh=blh;
         sesumh=sqrt(diag(vch));
         bsuml=bll;
         sesuml=sqrt(diag(vcl));
      endif;
      if trial>=2;
         ystats=ystats|(meanc(y)~stdc(y)~minc(y)~maxc(y));
         xstats=xstats~(meanc(x[.,2:cols(x)])|stdc(x[.,2:cols(x)])|
                minc(x[.,2:cols(x)])|maxc(x[.,2:cols(x)]));
         bsum=bsum~b;
         sesum=sesum~(sqrt(diag(vc)));
         bsuma=bsuma~ba;
         sesuma=sesuma~sqrt(diag(vca));
         bsumh=bsumh~blh;
         sesumh=sesumh~sqrt(diag(vch));
         bsuml=bsuml~bll;
         sesuml=sesuml~sqrt(diag(vcl));
      endif;

/* vectors of first difference values for the j^th independent var */
         bns=b[1:cols(x),.];                @ nxk matrix @
         hns=b[cols(x)+1:rows(b),.];   @ nxk matrix @

/* calculate model mse */

         ahat=exp(x*bns);
	 bhat=exp(x*hns);
         yhatns=ahat./(ahat+bhat);
	 bnsmse=sumc((y-yhatns)^2)/n;
	 
/* get predicted alpha and beta values with the j^th independent var 
   at +/- 1 sd and min/max values, while other x's held at their means */
         alns=exp(xl'*bns);                          @ nxk-1 matrix @
         blns=exp(xl'*hns);
         ahns=exp(xh'*bns);
         bhns=exp(xh'*hns);

/* get predicted values of mean and variance with the j^th independent var
   at +/- sd and min/max values, while other x's held at their means */
         meanlns=alns./(alns+blns);                    @ nxk-1 matrix @
         meanhns=ahns./(ahns+bhns);
         varlns=(alns.*blns)./(((alns+blns)^2).*(alns+blns+1));
         varhns=(ahns.*bhns)./(((ahns+bhns)^2).*(ahns+bhns+1));

/* produce first differences for the predicted mean and variance with 
   the j^th indep var at +/- sd, while other x's held at their means */
         fdmns=meanhns-meanlns;
         fdvns=varhns-varlns;

/* produce vector for average first differences */

         if trial==1 and ret==0; 
            totmns=fdmns;  @ k-1X1 matrix for mean @
            totvns=fdvns;  @ k-1X1 matrix for variance @
	    totbmse=bnsmse;
         endif;

         if trial>=2 and ret==0; 
            totmns=totmns~fdmns;  @ k-1X(trial) matrix for mean @
            totvns=totvns~fdvns;  @ k-1X(trial) matrix for variance @
	    totbmse=totbmse|bnsmse;
         endif;

/* vectors of first difference values for the j^th independent var */
         betasim=(rndmn(b,vc,1000));  @ nxk matrix, where n=# of draws @
         be=betasim[.,1:cols(x)];                @ nxk matrix @
         h=betasim[.,cols(x)+1:cols(betasim)];   @ nxk matrix @

/* get predicted alpha and beta values with the j^th independent var 
   at +/- 1 sd and min/max values, while other x's held at their means */
         al=exp(be*xl);                          @ nxk-1 matrix @
         bl=exp(h*xl);
         ah=exp(be*xh);
         bh=exp(h*xh);

/* get predicted values of mean and variance with the j^th independent var
   at +/- sd and min/max values, while other x's held at their means */
         meanlow=al./(al+bl);                    @ nxk-1 matrix @
         meanhigh=ah./(ah+bh);
         varlow=(al.*bl)./(((al+bl)^2).*(al+bl+1));
         varhigh=(ah.*bh)./(((ah+bh)^2).*(ah+bh+1));

/* produce first differences for the predicted mean and variance with 
   the j^th indep var at +/- sd, while other x's held at their means */
         fdm=meanhigh-meanlow;
         fdv=varhigh-varlow;

/* produce vector for average first differences */

         if trial==1 and ret==0; 
            totm=meanc(fdm);  @ k-1X1 matrix for mean @
            totv=meanc(fdv);  @ k-1X1 matrix for variance @
            totms=stdc(fdm);  @ k-1X1 matrix for mean @
            totvs=stdc(fdv);  @ k-1X1 matrix for variance @
         endif;

         if trial>=2 and ret==0; 
            totm=totm~meanc(fdm);  @ k-1X(trial) matrix for mean @
            totv=totv~meanc(fdv);  @ k-1X(trial) matrix for variance @
            totms=totms~stdc(fdm);  @ k-1X(trial) matrix for mean @
            totvs=totvs~stdc(fdv);  @ k-1X(trial) matrix for variance @
         endif;

/* get relevant stats for alternative method of beta estimation */
/* vectors of first difference values for the j^th independent var */
         bnsa=ba[1:cols(x),.];                @ nxk matrix @
         hnsa=ba[cols(x)+1:rows(ba),.];   @ nxk matrix @

/* calculate mse for alternative beta estimation */
  
         yhatnsa=exp(x*bnsa)./(1+exp(x*bnsa));
 	 bnsamse=sumc((y-yhatnsa)^2)/n;
	 
/* get predicted values of mean and variance with the j^th independent var
   at +/- sd and min/max values, while other x's held at their means */
         meanlnsa=(exp(xl'*bnsa)./(1+exp(xl'*bnsa)))';   @ nxk-1 matrix @
         meanhnsa=(exp(xh'*bnsa)./(1+exp(xh'*bnsa)))';
         dlsa=(1/(exp(xl'*hnsa)+1))';
         dhsa=(1/(exp(xh'*hnsa)+1))';
         varlnsa=meanlnsa.*(1-meanlnsa).*dlsa;
         varhnsa=meanhnsa.*(1-meanhnsa).*dhsa;

/* produce first differences for the predicted mean and variance with 
   the j^th indep var at +/- sd, while other x's held at their means */
         fdmnsa=meanhnsa-meanlnsa;
         fdvnsa=varhnsa-varlnsa;

/* produce vector for average first differences */

         if trial==1 and ret==0; 
            totmnsa=fdmnsa;  @ k-1X1 matrix for mean @
            totvnsa=fdvnsa;  @ k-1X1 matrix for variance @
	    totbamse=bnsamse;
         endif;

         if trial>=2 and ret==0; 
            totmnsa=totmnsa~fdmnsa;  @ k-1X(trial) matrix for mean @
            totvnsa=totvnsa~fdvnsa;  @ k-1X(trial) matrix for var @
	    totbamse=totbamse|bnsamse;
         endif;

/* vectors of first difference values for the j^th independent var */
         betasima=(rndmn(ba,vca,1000));  @ nxk matrix, where n=# of draws @
         bea=betasima[.,1:cols(x)];                @ nxk matrix @
         ha=betasima[.,cols(x)+1:cols(betasima)];   @ nxk matrix @

/* get predicted values of mean and variance with the j^th independent var
   at +/- sd and min/max values, while other x's held at their means */
         meanla=(exp(bea*xl)./(1+exp(bea*xl)));   @ nxk-1 matrix @
         meanha=(exp(bea*xh)./(1+exp(bea*xh)));
         dl=(1/(exp(ha*xl)+1));
         dh=(1/(exp(ha*xh)+1));
         varla=meanla.*(1-meanla).*dl;
         varha=meanha.*(1-meanha).*dh;

/* produce first differences for the predicted mean and variance with 
   the j^th indep var at +/- sd, while other x's held at their means */
         fdma=meanha-meanla;
         fdva=varha-varla;

/* produce vector for average first differences */

         if trial==1 and ret==0; 
            totma=meanc(fdma);  @ k-1X1 matrix for mean @
            totva=meanc(fdva);  @ k-1X1 matrix for variance @
            totmsa=stdc(fdma);
            totvsa=stdc(fdva);
         endif;

         if trial>=2 and ret==0; 
            totma=totma~meanc(fdma);  @ k-1X(trial) matrix for mean @
            totva=totva~meanc(fdva);  @ k-1X(trial) matrix for variance @
            totmsa=totmsa~stdc(fdma);
            totvsa=totvsa~stdc(fdva);
         endif;

/* set parameter vectors (linear heteroskeadastic) */
         bhh=blh[1:cols(x),1];
         hh=trimr(blh,rows(bhh),0);

/* calculate mse */

         yhatlh=x*bhh;
	 lhmse=sumc((y-yhatlh)^2)/n;
	 
/* get predicted means */
         elh=xl'*bhh;
         ehh=xh'*bhh;

/* get predicted variances */
         vlh=exp(xl'*hh);
         vhh=exp(xh'*hh);

/* produce first differences */
         if trial==1 and ret==0;
         totfdmh=ehh-elh;
         totfdvh=vhh-vlh;
	 totlhmse=lhmse;
         elseif trial>=2 and ret==0;
         totfdmh=totfdmh~(ehh-elh);
         totfdvh=totfdvh~(vhh-vlh);
	 totlhmse=totlhmse|lhmse;
         endif;

/* set parameter vectors (linear logit tranformation) */
         blo=bll[1:cols(x),1];
         hl=trimr(bll,rows(blo),0);
	 
/* calculate mse for linear het log trans */
  
         yhatlhlt=exp(x*blo)./(1+exp(x*blo));
	 lhltmse=sumc((y-yhatlhlt)^2)/n;

/* get predicted means */
         eol=exp(xl'*blo)./(1+exp(xl'*blo));
         eoh=exp(xh'*blo)./(1+exp(xh'*blo));

/* get predicted variances */
         vol=((exp(xl'*blo)./(1+exp(xl'*blo))^2)^2).*exp(xl'*hl);
         voh=((exp(xh'*blo)./(1+exp(xh'*blo))^2)^2).*exp(xh'*hl);

/* produce first differences */
         fdml=eoh-eol;
         fdvl=voh-vol;

/* produce first differences */
         if trial==1 and ret==0;
         totfdml=fdml;
         totfdvl=fdvl;
	 totthmse=lhltmse;
         elseif trial>=2 and ret==0;
         totfdml=totfdml~fdml;
         totfdvl=totfdvl~fdvl;
	 totthmse=totthmse|lhltmse;
         endif;

/* run ols */

{ vnam,m,bols,stb,vcols,stderr,sigma,cx,rsq,resid,dbw } = ols(0,y,x);

/* calculate mse */

      yhat=x*bols;
      olsmse=sumc((y-yhat)^2)/n;

/* obtain the overconfidence scores for the OLS parameters */

      if trial==1;
         blsum=bols;
         selsum=stderr;
      endif;
      if trial>=2;
         blsum=blsum~bols;
         selsum=selsum~stderr;
      endif;

/* obtain the first differences */

olsl=xl'*bols;
olsh=xh'*bols;
fdols=olsh-olsl;

      if trial==1;
         fdolsum=fdols;
	 totomse=olsmse;
      endif;
      if trial>=2;
         fdolsum=fdolsum~fdols;
	 totomse=totomse|olsmse;
      endif;

/* run ols (with logit trans) */

yl=ln(y./(1-y));

{ vnam,m,bolsl,stb,vcols,stderrl,sigma,cx,rsq,resid,dbw } = ols(0,yl,x);

/* calculate mse */

      yhatlt=exp(x*bolsl)./(1+exp(x*bolsl));
      ltmse=sumc((y-yhatlt)^2)/n;

/* obtain the overconfidence scores for the OLS parameters */

      if trial==1;
         blsumo=bolsl;
         selsumo=stderrl;
      endif;
      if trial>=2;
         blsumo=blsumo~bolsl;
         selsumo=selsumo~stderrl;
      endif;

/* obtain the first differences */

olslo=exp(xl'*bolsl)./(1+exp(xl'*bolsl));
olsho=exp(xh'*bolsl)./(1+exp(xh'*bolsl));
fdolsl=olsho-olslo;

      if trial==1;
         fdolsuml=fdolsl;
	 tototmse=ltmse;
      endif;
      if trial>=2;
         fdolsuml=fdolsuml~fdolsl;
	 tototmse=tototmse|ltmse;
      endif;

/* run ols (with log trans) */

ly=ln(y+.01);

{ vnam,m,bolsll,stb,vcols,stderrll,sigma,cx,rsq,resid,dbw } = ols(0,ly,x);

/* calculate mse */

      yhatll=exp(x*bolsll)-0.01;
      ltmsel=sumc((y-yhatll)^2)/n;

/* obtain the overconfidence scores for the OLS parameters */

      if trial==1;
         blsumol=bolsll;
         selsumol=stderrll;
      endif;
      if trial>=2;
         blsumol=blsumol~bolsll;
         selsumol=selsumol~stderrll;
      endif;

/* obtain the first differences */

olslol=exp(xl'*bolsll)-0.01;
olshol=exp(xh'*bolsll)-0.01;
fdolsll=olshol-olslol;

      if trial==1;
         fdsumll=fdolsll;
	 totmsell=ltmsel;
      endif;
      if trial>=2;
         fdsumll=fdsumll~fdolsll;
	 totmsell=totmsell|ltmsel;
      endif;

/* reset fails and add a trial for a successful trial */
      fails=0;
      trial=trial+1; 
      clear y;

   endif;

endo;

/* check that parameters are accurate */
print "descriptive stats for standard beta parameters";
meanc(bsum')~stdc(bsum'); 
print " ";
/* check the mse for the std beta */
  
print "standard beta mse=";;meanc(totbmse);
print " ";
print "alternative beta mse=";;meanc(totbamse);
print " ";
print "ols mse=";;meanc(totomse);
print " ";
print "logit trans ols mse=";;meanc(tototmse);
print " ";
print "het ols mse=";;meanc(totlhmse);
print " ";
print "logit trans het ols mse=";;meanc(totthmse);
print " ";
print "log trans het ols mse=";;meanc(totmsell);
print " ";

/* check mean and std. dev of first differences */

print"first differences for mean and variance for std beta (sim)";
meanc(totm')~stdc(totm');
meanc(totv')~stdc(totv');
print " ";

print "first differences for mean and variance for std beta (nosim)";
meanc(totmns')~stdc(totmns');
meanc(totvns')~stdc(totvns');
print " ";

print "first differences for mean and variance for alt beta (sim)";
meanc(totma')~stdc(totma');
meanc(totva')~stdc(totva');
print " ";

print "first differences for mean and variance for alt beta (nosim)";
meanc(totmnsa')~stdc(totmnsa');
meanc(totvnsa')~stdc(totvnsa');
print " ";

/* check for mean and std. dev. of first difference for OLS */
print"first differences for mean ols";
meanc(fdolsum')~stdc(fdolsum');
print " ";

/* check for mean and std. dev. of first difference for OLS (log trans) */
print"first differences for mean ols (logit trans)";
meanc(fdolsuml')~stdc(fdolsuml');
print " ";

print "first differences for mean and variance for het normal";
meanc(totfdmh')~stdc(totfdmh');
meanc(totfdvh')~stdc(totfdvh');
print " ";

print "first differences for mean and variance for logit trans.";
meanc(totfdml')~stdc(totfdml');
meanc(totfdvl')~stdc(totfdvl');
print " ";

/* check for mean and std. dev. of first difference for OLS (log trans) */
print"first differences for mean ols (log trans)";
meanc(fdsumll')~stdc(fdsumll');
print " ";

/* efficiency of expected values */
/* compare efficiency of standard beta with alternative beta */

effnba=((totmnsa-fdmt)^2);
effdba=((totmns-fdmt)^2);
effalt=100*(sqrt(sumc(effnba')))./(sqrt(sumc(effdba')));
print "efficiency of standard beta vs. alternative beta=";;effalt;
print " ";

/* compare efficiency of standard beta with ols */

effno=((fdolsum-fdmt)^2);
effdb=((totmns-fdmt)^2);
effols=100*(sqrt(sumc(effno')))./(sqrt(sumc(effdb')));
print "efficiency of beta vs. ols=";;effols;
print " ";

/* compare efficiency of standard beta with ols */

effnl=((fdolsuml-fdmt)^2);
effdb=((totmns-fdmt)^2);
effolsl=100*(sqrt(sumc(effnl')))./(sqrt(sumc(effdb')));
print "efficiency of beta vs. ols (log trans)=";;effolsl;
print " ";

/* compare efficiency of standard beta with heteroskedasic normal */

effnbh=((totfdmh-fdmt)^2);
effdba=((totmns-fdmt)^2);
effhet=100*(sqrt(sumc(effnbh')))./(sqrt(sumc(effdba')));
print "efficiency of standard beta vs. het normal=";;effhet;
print " ";

/* compare efficiency of standard beta with heteroskedasic normal */

effnbll=((totfdml-fdmt)^2);
effdbal=((totmns-fdmt)^2);
efflog=100*(sqrt(sumc(effnbll')))./(sqrt(sumc(effdbal')));
print "efficiency of standard beta vs. het normal (log trans)=";;efflog;
print " ";

/* compare efficiency of standard beta with log trans normal */

effnbl=((fdsumll-fdmt)^2);
effdba=((totmns-fdmt)^2);
efflogl=100*(sqrt(sumc(effnbl')))./(sqrt(sumc(effdba')));
print "efficiency of standard beta vs. OLS (log trans)=";;efflogl;
print " ";

/* efficiency for variance term */
/* compare efficiency of standard beta with alternative beta */

veffnba=((totvnsa-fdvt)^2);
veffdba=((totvns-fdvt)^2);
veffalt=100*(sqrt(sumc(veffnba')))./(sqrt(sumc(veffdba')));
print "efficiency of standard beta vs. alternative beta=";;veffalt;
print " ";

/* compare efficiency of standard beta with heteroskedasic normal */

veffnbh=((totfdvh-fdvt)^2);
veffdba=((totvns-fdvt)^2);
veffhet=100*(sqrt(sumc(veffnbh')))./(sqrt(sumc(veffdba')));
print "efficiency of standard beta vs. het normal=";;veffhet;
print " ";

/* compare efficiency of standard beta with heteroskedasic normal */

veffnbll=((totfdvl-fdvt)^2);
veffdbal=((totvns-fdvt)^2);
vefflog=100*(sqrt(sumc(veffnbll')))./(sqrt(sumc(veffdbal')));
print "efficiency of standard beta vs. het normal (log trans)=";;vefflog;
print " ";

/* check standard beta parameters' overconfidence */

npoca=(bsum-meanc(bsum'))^2;
npoc=sqrt(sumc(npoca'));
dpoca=sesum^2;
dpoc=sqrt(sumc(dpoca'));
ocp=100*(npoc./dpoc);
print "standard beta parms' overconfidence=";;ocp;
clear npoca,npoc,dpoca,dpoc;
print " ";

/* check alternative beta parameters' overconfidence */

npoca=(bsuma-meanc(bsuma'))^2;
npoc=sqrt(sumc(npoca'));
dpoca=sesuma^2;
dpoc=sqrt(sumc(dpoca'));
ocpa=100*(npoc./dpoc);
print "alternative beta parms' overconfidence=";;ocpa;
clear npoca,npoc,dpoca,dpoc;
print " ";

/* check ols parameters' overconfidence */

npoca=(blsum-meanc(blsum'))^2;
npoc=sqrt(sumc(npoca'));
dpoca=selsum^2;
dpoc=sqrt(sumc(dpoca'));
oco=100*(npoc./dpoc);
print "ols parms' overconfidence=";;oco;
clear npoca,npoc,dpoca,dpoc;
print " ";

/* check ols (log. trans) parameters' overconfidence */

npoca=(blsumo-meanc(blsumo'))^2;
npoc=sqrt(sumc(npoca'));
dpoca=selsumo^2;
dpoc=sqrt(sumc(dpoca'));
ocol=100*(npoc./dpoc);
print "ols (log. trans) parms' overconfidence=";;ocol;
clear npoca,npoc,dpoca,dpoc;
print " ";

/* check ols (log trans) parameters' overconfidence */

npoca=(blsumol-meanc(blsumol'))^2;
npoc=sqrt(sumc(npoca'));
dpoca=selsumol^2;
dpoc=sqrt(sumc(dpoca'));
ocoll=100*(npoc./dpoc);
print "ols (log trans) parms' overconfidence=";;ocoll;
clear npoca,npoc,dpoca,dpoc;
print " ";

/* check heteroskedastic normal parameters' overconfidence */

npoca=(bsumh-meanc(bsumh'))^2;
npoc=sqrt(sumc(npoca'));
dpoca=sesumh^2;
dpoc=sqrt(sumc(dpoca'));
ocph=100*(npoc./dpoc);
print "heteroskedastic normal parms' overconfidence=";;ocph;
clear npoca,npoc,dpoca,dpoc;
print " ";

/* check heteroskedastic logit transform parameters' overconfidence */

npoca=(bsuml-meanc(bsuml'))^2;
npoc=sqrt(sumc(npoca'));
dpoca=sesuml^2;
dpoc=sqrt(sumc(dpoca'));
ocph=100*(npoc./dpoc);
print "heteroskedastic log. trans. parms' overconfidence=";;ocph;
clear npoca,npoc,dpoca,dpoc;
print " ";

/* check overconfidence for the first differences (mean) */

print "standard beta's mean first difference overconfidence";
nmoca=(totm-meanc(totm'))^2;
nmoc=sqrt(sumc(nmoca'));
dmoca=totms^2;
dmoc=sqrt(sumc(dmoca'));
ocm=100*(nmoc./dmoc);
ocm;
print " ";

/* check overconfidence for the first differences (variance) */

print "standard beta's variance first difference overconfidence";
nvoca=(totv-meanc(totv'))^2;
nvoc=sqrt(sumc(nvoca'));
dvoca=totvs^2;
dvoc=sqrt(sumc(dvoca'));
ocv=100*(nvoc./dvoc);
ocv;
print " ";

/* check overconfidence for the first differences (mean) */

print "alternative beta's mean first difference overconfidence";
nmocaa=(totma-meanc(totma'))^2;
nmoca=sqrt(sumc(nmocaa'));
dmocaa=totmsa^2;
dmoca=sqrt(sumc(dmocaa'));
ocma=100*(nmoca./dmoca);
ocma;
print " ";

/* check overconfidence for the first differences (variance) */

print "alternative beta's variance first difference overconfidence";
nvocaa=(totva-meanc(totva'))^2;
nvoca=sqrt(sumc(nvocaa'));
dvocaa=totvsa^2;
dvoca=sqrt(sumc(dvocaa'));
ocva=100*(nvoca./dvoca);
ocva;
print " ";

/* check percentage of parameters within 95% conf int. */

print "check standard beta's percentage within 95% conf int";
cib=abs(bsum-meanc(bsum')).<(1.96*sesum);
sumc(cib')/(trial-1);
print " ";

print "check alternative beta's percentage within 95% conf int";
ciba=abs(bsuma-meanc(bsuma')).<(1.96*sesuma);
sumc(ciba')/(trial-1);
print " ";

print "check ols percentage within 95% conf int";
cio=abs(blsum-meanc(blsum')).<(1.96*selsum);
sumc(cio')/(trial-1);
print " ";

print "check logit trans/ols percentage within 95% conf int";
cit=abs(blsumo-meanc(blsumo')).<(1.96*selsumo);
sumc(cit')/(trial-1);
print " ";

print "check het ols percentage within 95% conf int";
cih=abs(bsumh-meanc(bsumh')).<(1.96*sesumh);
sumc(cih')/(trial-1);
print " ";

print "check logit trans/ols percentage within 95% conf int";
cil=abs(bsuml-meanc(bsuml')).<(1.96*sesuml);
sumc(cil')/(trial-1);
print " ";

print "check logit trans/ols percentage within 95% conf int";
cill=abs(blsumol-meanc(blsumol')).<(1.96*selsumol);
sumc(cill')/(trial-1);
print " ";

/* get information about percentage of significant parameters */
/* mean values */
print "check percentage of sig parms (standard beta)";
tstatbm=(abs(totm./totms)).>1.96;
tstatbv=(abs(totv./totvs)).>1.96;
sumc(tstatbm')/(trial-1);
sumc(tstatbv')/(trial-1);
print " ";

print "check percentage of sig parms (alternate beta)";
tstatba=(abs(bsuma./sesuma)).>1.96;
sumc(tstatba')/(trial-1);
print " ";

print "check percentage of sig parms (basic ols)";
tstatbo=(abs(blsum./selsum)).>1.96;
sumc(tstatbo')/(trial-1);
print " ";

print "check percentage of sig parms (logit trans ols)";
tstatbl=(abs(blsumo./selsumo)).>1.96;
sumc(tstatbl')/(trial-1);
print " ";

print "check percentage of sig parms (het ols)";
tstatboh=(abs(bsumh./sesumh)).>1.96;
sumc(tstatboh')/(trial-1);
print " ";

print "check percentage of sig parms (het ols - logit trans)";
tstatblh=(abs(bsuml./sesuml)).>1.96;
sumc(tstatblh')/(trial-1);
print " ";

print "check percentage of sig parms (ols -- log tranform)";
tstatbm=(abs(blsumol./selsumol)).>1.96;
sumc(tstatbm')/(trial-1);
print " ";


/* obtain information about y's */
print "descriptive stats for y's for successful runs";
meanc(ystats);
print "descriptive stats for y's for failed runs";
meanc(ystatf);

/* obtain information about x's */
print "descriptive stats for x's for successful runs";
meanc(xstats');
print "descriptive stats for x's for failed runs";
meanc(xstatf');

print "number of failed runs";
crash;
output off;